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ABSTRACT 


A  major  contributor  to  the  error  in  regional  and  teleseismic  event  locations  is  travel-time  prediction  error, 
attributable  to  velocity  anomalies  in  the  real  Earth  that  are  not  represented  in  the  reference  Earth  model  an  event 
locator  uses  for  travel-time  calculation.  The  prediction  errors  at  two  stations  will,  in  general,  be  correlated  depending 
on  the  spatial  coherency  of  the  velocity  anomalies  compared  to  the  distance  between  the  stations.  The  motivation  for 
this  project  is  the  observation  by  previous  workers  that  large  event  location  errors  can  be  induced  when  prediction 
error  correlations  are  ignored  and  the  observing  network  of  stations  is  far  from  uniform.  Recently,  event  location 
algorithms  have  been  generalized  to  accept  a  non-diagonal  covariance  matrix  for  data  errors  as  a  mechanism  for 
accommodating  correlated  errors  in  travel-time  predictions.  Our  project  is  addressing  what  to  input  for  the 
covariance  matrix,  based  on  considerations  of  seismic  wave  propagation  and  the  statistical  characterization  of  the 
Earth’s  velocity  heterogeneity.  Specifically,  we  are  developing  numerical  algorithms  that  generate  a  travel-time 
covariance  matrix  for  a  network  of  stations,  as  a  function  of  the  event  location,  by  integrating  a  velocity  covariance 
function,  as  defined  in  geostatistical  analysis,  with  the  travel-time  sensitivity  kernels  appropriate  for  the 
event-station  paths.  This  paper  describes  two  alternative  techniques  for  performing  this  calculation,  which  are 
compared  in  a  numerical  example.  A  second  example  of  our  approach  investigates  the  variation  of  travel-time 
variances  and  correlations  over  regional  and  near-teleseismic  distances,  computed  using  ray  sensitivities  for  a 
ID  reference  Earth  model.  This  example  demonstrates  a  strong  dependence  on  ray  parameter  that  agrees 
qualitatively  with  empirical  observations. 
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OBJECTIVE 


To  achieve  accurate  event  locations  and  estimates  of  location  uncertainty,  seismic  event  location  algorithms  must 
account  for  two  types  of  errors:  pick  errors  in  the  measured  arrival  times  of  the  seismic  phases  observed  at  various 
stations,  and  model  errors  incurred  in  predicting  travel  time  for  the  observed  station/phase  combinations  as  a 
function  of  the  event  location.  Model  errors  are  attributable  to  velocity  anomalies  in  the  Earth  that  are  not  rendered 
in  the  velocity  model  used  for  travel-time  prediction.  When  the  observing  stations  are  sufficiently  close  to  one 
another,  their  model  errors  are  expected  to  be  correlated.  While  it  is  now  standard  practice  for  event  location 
algorithms  used  in  nuclear  monitoring  to  assign  error  variances  that  include  the  contribution  from  model  prediction 
errors,  the  algorithms  generally  set  the  error  covariances  to  zero;  correlations  are  ignored.  Doing  so  can  seriously 
degrade  location  accuracy  when  the  distribution  of  seismic  stations  is  far  from  uniform  (e.g.  Chang  et  al.,  1983; 

Yang  et  al.  2004). 

This  project  is  investigating  the  phenomenon  of  correlated  travel-time  prediction  errors  from  a  physical  point  of 
view,  making  use  of  the  principles  of  wave  propagation  through  heterogeneous  media  and  our  knowledge  of  the 
statistical  properties  of  the  Earth’s  heterogeneity.  Specifically,  we  are  developing  techniques  to  calculate  covariance 
matrices  of  travel-time  model  errors  by  integrating  travel-time  sensitivity  kernels  with  plausible  correlation 
functions  describing  the  Earth's  velocity  heterogeneity.  Our  preliminary  efforts  in  this  development  were  reported  by 
Rodi  and  Myers  (2006).  Here  we  summarize  our  approach  and  report  on  our  recent  efforts. 

RESEARCH  ACCOMPITSHED 


Formulation  of  Approach 

Given  a  set  of  n  arrival  time  data  from  an  event,  one  can  decompose  the  data  errors  as  (for  ;  =  !,...,«) 

(1) 

where  the  first  term  is  the  pick,  or  measurement,  error  and  the  second  term  is  the  model  error,  or  error  in  the  travel- 
time  prediction  from  a  reference  Earth  velocity  model.  Typically,  event  location  algorithms  assume  that  the 
variances  assigned  to  the  arrival  time  data,  include  a  contribution  for  each  type  of  error: 


2  2  2 

0-,  =(7p,.  +  (7„,.. 


(2) 


However,  the  errors  for  different  i  are  assumed  uncorrelated.  We  are  addressing  the  problem  of  calculating  a  full 
covariance  matrix  for  model  errors,  having  components  defined  by 

(3) 


E[  ]  denotes  the  expectation  operator. 

To  explain  our  approach  to  covariance  modeling,  we  consider  only  P-wave  arrivals.  Let  Vo(x)  denote  the  P  velocity 
function  for  the  reference  model,  and  V£(x)  denote  the  Earth’s  true  (and  unknown)  velocity  function.  Then,  model 
errors  can  be  linked  to  the  slowness  difference,  which  we  denote  m(x): 

«(x)  =  V£  (x)-Vo‘(x).  (4) 

We  assume  that  the  travel-time  dependence  on  slowness  can  be  adequately  approximated  as  linear,  allowing  us  to 
express  the  model  error  in  the  fth  datum  as 


em,i=\  ai{^)m{^)d^,  (5) 

where  a,(x)  is  the  first-order  sensitivity  kernel  of  the  /th  travel-time  functional,  as  evaluated  at  vq.  In  the 
high-frequency  limit,  this  kernel  is  concentrated  along  the  geometrical  ray  connecting  the  event  and  station 
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locations.  For  finite  frequency,  it  is  spatially  distributed  around  this  ray.  In  either  case,  we  point  out  that  the  model 
error  is  a  function  of  the  event  and  station  locations. 


We  consider  w(x)  to  be  a  Gaussian  random  field,  having  zero  mean  and  a  specified  covariance  between  any  two 
points,  described  by  a  covariance  operator  C(x,x') : 

(6) 
(7) 


E[»i(x)]=  0 

E  [ot  (x  )m  (x')]  =  C  (x ,  x'). 


The  covariance  between  two  model  errors  is  then  given  by 

o-m.y  =  J  t/xa,.(x)J  Jx'C(x,x')a,(x'). 


(8) 


Equation  (8)  implies  that  the  extent  to  which  two  model  errors  are  correlated  depends  on  the  spatial  relationship 
between  their  sensitivity  kernels  in  relation  to  the  correlation  structure  of  the  velocity  field. 


Velocity  Model  Covariance  Operators 

When  m(x)  is  a  stationary  random  field,  C(x,x')  depends  only  on  the  difference  between  the  points  x  and  x'.  For 
example,  if  m(x)  is  isotropic  with  a  correlation  function  of  the  exponential  type  (Deutsch  and  Joumel,  1998),  then 


C(x,x')  =  a^  exn  - 


A 


(9) 


where  is  the  variance  of  the  slowness  field  and  A  is  a  correlation  length.  Flowever,  the  Earth’s  velocity 

heterogeneity  is  not  well-approximated  as  stationary  over  the  spatial  scales  affecting  regional  and  teleseismic  waves. 
It  is  difficult  to  generalize  an  operator  such  as  the  exponential  form  above  to  allow  for  nonstationarity  while 
preserving  the  essential  properties  of  C  (x,  x') ,  in  particular  its  being  a  positive-definite  operator. 

A  flexible  approach  to  characterizing  nonstationary  random  fields  is  to  specify  the  covariance  operator  indirectly 
through  its  operator  inverse,  which  we  denote  D,  such  that 

Z)C(x,x')  =  <5(x-x').  (10) 

If  we  take  Z)  to  be  a  differential  operator,  then  C(x,x')  is  its  Green's  function.  Rodi  et  al.  (2003)  implemented  this 
approach  with  D  as 


D  = 


const 


<J(x)- 


(2/ -3) 


A  V' 

\r"  ' 


dz\ 


I  >2. 


(11) 


Flere,  /  is  an  integer  specifying  the  order  of  the  operator;  Vj*  is  the  horizontal  Laplacian  operator;  Ai  and  A2  are 
horizontal  and  vertical  correlation  lengths,  respectively;  and  is  the  variance.  The  exponential-type  covariance 
operator  of  equation  (9)  is  given  by  /  =  2,  ignoring  boundaries. 

In  previous  work  (Rodi  and  Myers,  2006)  we  have  implemented  our  numerical  algorithms  for  covariance  modeling 
with  D  being  nonstationary  only  by  virtue  of  the  boundary  at  the  Earth’s  surface  and  velocity  discontinuities 
(e.g.,  the  Moho),  across  which  C(x,x')  is  allowed  to  “de-correlate”,  i.e.  C(x,x')  =  0  when  x  and  x'  are  on  opposing 
sides  of  the  discontinuity.  One  of  the  accomplishments  of  the  current  year  has  been  to  implement  a  more  general 
form  of  nonstationarity  whereby  crand  A2  (vertical  correlation  distance)  are  allowed  to  vary  with  depth. 
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Covariance  Modeling  Algorithms 

The  numerical  techniques  we  have  developed  apply  to  a  discrete  model  parameterization.  We  thus  replace  the  model 
function  m(x)  with  a  A:  X 1  vector  m;  the  sensitivity  kernels  a,(x)  with  A:  X 1  vectors  a,;  and  the  operators  C  and  D 
with  A:  X  A:  matrices  C  and  D.  We  require  that 


DC  =  I, 


(12) 


where  I  is  the  kxk  identity  matrix. 

In  the  discrete  formulation,  the  double  integral  of  equation  (8)  becomes  the  quadratic  matrix  expression 


(13) 


It  is  also  useful  to  express  the  full  n  x  «  model-error  covariance  matrix,  which  we  denote  Vm  and  whose  elements 
are  the  .  To  do  this,  we  first  define  the  nxk  sensitivity  matrix  having  the  sensitivity  vectors  as  its  rows: 

A''  =  (a,  a,  •••  a„).  (14) 


We  can  then  write 


V„  =  ACA^  (15) 

We  now  outline  two  numerical  methods  we  have  developed  for  computing  Both  assume  that  D  is  specified 
directly,  with  C  then  being  given  implicitly  as  the  inverse  of  D.  Both  methods  require  that  the  sensitivity  vectors,  a„ 
be  given  explicitly. 

Explicit  Method 

This  method  entails  solving  a  linear  system  for  each  j,  given  by 

DUy  =  ay. 

Given  the  solution  vector  Uy,  the  elements  of  one  column  of  the  model-error  covariance  matrix  can  be  found  as 

<T  .  =  a'^M 1  = !,...,« 

since,  analytically,  U/=  Cay.  Repeating  for  each j=  the  full  matrix  Vm  is  generated. 

Being  an  approximation  to  a  differential  operator,  the  matrix  D  is  highly  sparse.  Therefore,  our  implementation  of 
this  method  solves  (16)  with  a  conjugate  gradients  technique. 

Theoretically,  the  matrix  obtained  with  the  explicit  method  will  be  symmetric  since  ufa  =  u^a .  if  the  solution 
vectors  Uy  are  exact.  Numerically,  however,  this  will  not  be  the  case  exactly.  Therefore,  we  average  the  two 
numerical  estimates  of  (Tmy  (I  j)  to  generate  a  symmetric  approximation  to  Vm. 

Implicit  Method 

This  method  was  described  by  Rodi  and  Myers  (2006)  and  we  summarize  it  here  for  completeness. 

For  each  j  we  solve  a  minimization  problem  for  a  vector  in  the  model  space,  Wy: 

(n^-Aw^)  Vp‘(n, -Aw  .^-t-w^Dw^  =  minimum  (18) 


(16) 

(17) 
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where  Vp  is  some  assumed,  diagonal  eovarianee  matrix  for  piek  errors,  and  n,  is  the yth  eolumn  of  the  n  x  n  identity 
matrix.  Again,  we  eompute  the  solution  W/  with  a  eonjugate  gradients  teehnique.  Given  the  solution,  we  form  the 
residual  veetor 


r .  =  n ,  -  Aw , . 

/  /  / 

Repeating  this  proeess  for  eaeh y  =  1,  . . .,  n,  we  generate  a  residual  matrix  eolumn-by-eolumn: 

R=(ri  G  •••  r„). 


Rodi  and  Myers  (2006)  showed  that 


R=V^(ACA^  +  F^y‘. 


(19) 


(20) 


(21) 


Therefore,  two  useful  quantities  can  be  generated  from  R.  The  first  is  the  inverse  of  the  total  covariance  matrix: 


(v„  +  F,y  =  F;'R 


(22) 


This  matrix  can  be  used  by  an  event  location  algorithm  in  the  calculation  of  a  misfit  function,  assuming  Vp  is 
appropriate.  Second,  we  have  the  object  of  this  project,  the  model-error  covariance  matrix: 

v„  =  (v;'Ry‘-v^.  (23) 

As  before,  the  numerically  generated  obtained  from  this  formula  is  not  guaranteed  to  be  symmetric,  but  we  can 
average  it  with  its  transpose  to  make  it  so. 

For  the  purpose  of  calculating  by  the  implicit  method,  the  matrix  Vp  is  arbitrary.  We  have  found  that  setting  it  to 
a  fraction  of  the  identity  yields  accurate  results.  If  Vp  is  sufficiently  small,  we  can  ignore  it  and  directly  obtain  a 
third  quantity  which  can  be  useful  in  a  location  algorithm: 

V„‘  »  V^^'R.  (24) 


Comparison  of  Explicit  and  Implicit  Methods 

Rodi  and  Myers  (2006)  presented  an  example  of  a  travel-time  (model-error)  covariance  matrix  for  a  network  of 
12  regional  stations  observing  a  shallow  earthquake  in  northeastern  Iran.  Figure  1  displays  the  station  geometry  for 
this  example.  The  event  comes  from  the  EFIB  (Engdahl  et  al.,  1998)  data  base  and  is  one  of  the  events  used  by 
Reiter  and  Rodi  (2006)  in  their  body-wave  tomography  study  of  southern  Asia.  The  covariance  calculations  were 
performed  using  sensitivity  vectors  (discretized  kernels)  generated  by  Reiter  and  Rodi  (2006)  for  their  initial 
3D  velocity  model,  calculated  with  a  finite-difference  raytracing  technique.  The  results  we  show  here  have  been 
updated  using  a  more  recent  tomographic  model  (Rodi  and  Reiter,  2007)  and  revised  values  of  geostatistical 
parameters  for  the  Earth’s  velocity.  Velocity  variations  were  assigned  a  standard  deviation  (o)  of  3%  in  the  cmst  and 
1.5%  in  the  upper  mantle.  The  vertical  correlation  distance  (/I2)  in  the  crust  was  set  to  half  the  crustal  thickness, 
while  /I2  =  60  km  was  used  throughout  the  upper  mantle.  The  horizontal  correlation  distance  was  300  km  in  the  crust 
and  upper  mantle. 

We  display  two  aspects  of  the  12  x  12  model-error  covariance  matrices  (V„)  resulting  from  our  calculations.  Figure 
2  displays  the  12  standard  deviations  (square-root  of  the  diagonal  elements)  as  a  function  of  epicenfral  distance.  The 
66  independent  off-diagonal  elements,  normalized  as  correlation  coefficients,  are  plotted  as  a  function  of 
inter-station  azimuth  in  Figure  3.  The  distance  and  azimuth  dependence  of  the  covariances  behave  as  expected, 
qualitatively  at  least  validating  our  modeling  approach  (see  Rodi  and  Myers,  2006,  for  a  more  detailed  discussion  of 
the  results).  For  our  purposes  here,  we  point  out  that  the  two  numerical  techniques  we  applied,  which  are  compared 
in  each  figure,  produce  very  similar  results,  although  some  small  difference  can  be  seen. 
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Fignre  1.  Event/station  geometry  nsed  for  testing  onr  travel-time  covariance  modeling  algorithms.  The  event 
epicenter  is  marked  as  the  red  circle.  The  event  had  Pn  arrivals  for  12  stations,  marked  as  hlne 
circles. 
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Fignre  2.  Standard  deviation  of  model  error  as  a  fnnction  of  epicentral  distance  compnted  with  the  explicit 
method  (left)  and  implicit  method  (right).  The  event-station  geometry  is  shown  in  Fignre  1. 

Distance  Dependence  of  Travel-Time  Covariance 

Our  second  example  uses  a  linear  array  of  70  stations  equally  spaced  from  0.5°  to  35°  in  epicentral  distance  from  a 
nominal  event  location  at  a  common  azimuth.  The  resulting  70  x  70  covariance  matrix  then  reveals  the  dependence 
of  travel-time  variance  on  distance  and  of  travel-time  correlation  on  distance  and  distance  separation  between 
stations.  The  calculations  were  done  using  travel-time  sensitivities  for  the  akl35  ID  Earth  model  (Kennett  et  al., 
1995),  obtained  with  analytical  solutions  for  geometrical  rays  (e.g.,  Buland  and  Chapman,  1983)  and  transformed  to 
discrete  sensitivities  for  a  3D  model  parameterization.  The  covariance  matrix  computations  were  perfonued  with  our 
explicit  method. 

We  performed  the  calculations  under  two  assumptions  about  the  geostatistical  parameters  of  velocity  heterogeneity. 
In  each,  the  correlation  distances  were  the  same  as  in  the  previous  example  =  300  km,  =  17.5  km  in  the  crust 
and  60  km  in  the  mantle).  The  standard  deviation  of  velocity  variations  in  the  crust  was  also  common  to  both  cases: 
<T=  3%  as  before.  The  two  geostatistical  models  differed  in  the  velocity  standard  deviation  assigned  to  the  mantle. 
The  first  case  used  <7=1.5%  throughout  the  mantle.  The  second  case  had  <7=  2%  in  the  upper  mantle  (above  the 
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410-km  discontinuity)  and  (7=  1%  below  410  km.  These  parameters  were  chosen  to  illustrate  the  parameter 
dependence  and  are  not  necessarily  indicative  of  the  actual  geostatistical  properties  of  the  Earth. 


Figure  3.  Correlation  coefficient  between  model-errors  as  a  function  of  azimuth  difference  between  stations, 
computed  with  the  explicit  (left)  and  implicit  (right)  methods.  The  event-station  geometry  is  shown 
in  Figure  1. 

Figure  4  shows  the  standard  deviation  of  travel-time  model  errors  for  the  two  geostatistical  models,  each  plotted 
versus  epicentral  distance.  Distance-dependent  discontinuities  in  the  travel-time  standard  deviation  are  evident  at  the 
same  distances  in  both  models.  These  discontinuities  are  directly  attributed  to  discontinuities  in  ray  parameter, 
which  are  controlled  by  the  akl35  model.  From  approximately  2°  to  15°  travel-time  error  increases  consistently,  but 
not  strictly  linearly.  The  rate  of  increase  in  the  error  diminishes  with  distances,  as  rays  begin  to  average  over  several 
correlation  lengths  of  velocity  variations.  At  approximately  15°,  rays  begin  to  dive  deeper  into  the  upper  mantle  and 
travel  more  vertically,  resulting  in  more  averaging  over  the  shorter  correlation-length  anomalies  in  the  vertical 
dimension.  Several  breaks  in  the  error  structure  are  evident  as  rays  dive  below  the  410-km  velocity  discontinuity  and 
then  the  660-km  discontinuity.  Error  structure  stabilizes  at  teleseismic  distances  (>  24°)  when  rays  bottom  in  the 
lower  mantle  and  the  ray  parameter  remains  continuous.  Travel-time  error  is  lower  at  teleseismic  distance  because 
vertically  traveling  rays  in  the  more  strongly  heterogeneous  upper  mantle  average  over  more  wavelengths  of 
geologic  heterogeneity,  recalling  that  we  assumed  the  vertical  scale  of  heterogeneity  to  be  much  smaller  than  the 
horizontal  scale. 

The  travel-time  error  structure  shown  in  Figure  4  is  similar  to  that  derived  empirically  (e.g.,  Flanagan  et  al.,  2007). 
The  similarity  between  empirical  and  our  model-based  predictions  of  travel-time  error  as  a  function  of  distance  is 
encouraging.  Further,  these  results  provide  insight  into  the  non-intuitive  but  robust  observation  that  teleseismic 
travel-time  errors  with  long  ray  paths  are  smaller  than  regional  travel-time  errors  with  shorter  ray  paths. 

Figure  5  shows  the  correlation  matrix  corresponding  to  the  standard  deviations  shown  in  Figure  4.  Analogous  to  the 
breaks  in  travel-time  standard  deviation.  Figure  5  shows  a  roughly  block-diagonal  structure,  with  the  blocks 
delimited  by  the  same  crossover  distances  and  associated  jumps  in  ray  parameter.  These  results  demonstrate  that 
simple  estimates  of  travel-time  correlation  based  on  the  distance  between  stations  do  not  adequately  account  for 
correlations  between  travel-times  near  crossover  distances.  If  two  stations  straddle  a  break  in  ray  parameter,  then 
travel-time  residuals  are  likely  to  be  uncorrelated.  This  result  has  profound  implications  for  both  location  algorithms 
and  empirical  methods  (e.g.,  kriging)  that  make  use  of  the  travel-time  residual  covariance  matrix. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


Our  efforts  to  date  eonfirm  the  feasibility  of  using  a  model-based  approaeh  to  eompute  travel-time  residual 
eovarianees.  The  numerieal  teehniques  we  have  developed  are  quite  flexible  in  the  station  geometries  and  types  of 


Figure  4.  Standard  deviation  of  travel-time  model  error  as  a  function  of  epicentral  distance  computed  with 
different  geostatistical  parameters  for  velocity  heterogeneity.  Both  have  a  3%  velocity  standard 
deviation  (c^  in  the  crust.  The  results  on  the  left  correspond  to  <T=  1.5%  everywhere  in  the  mantle, 
while  those  on  the  right  assign  C7=  2%  in  the  upper  mantle  above  a  depth  of  410  km,  and  tT=  1% 
helow  410  km.  The  same  correlation  lengths  of  heterogeneity  were  used  in  both  cases  (see  text). 


Distance  (deg)  ) 


Figure  5.  Correlation  coefficient  between  model  errors  at  different  epicentral  distances  and  a  common 
azimuth,  computed  with  the  same  two  geostatistical  models  compared  in  Figure  4;  (7=  1.5% 
throughout  the  mantle  (left)  and  cr=  2%  above  410  km  and  1%  below  410  km  (right). 


Earth  models  they  ean  aeeommodate  (ID  and  3D),  providing  a  useful  tool  for  the  theoretieal  investigation  of  how 
travel-time  eovarianees  behave  in  a  variety  of  situations.  Applying  the  teehniques  under  even  simple  assumptions 
about  the  statistieal  properties  of  veloeity  anomalies  in  the  Earth  has  already  yielded  important  insights  into  previous 
empirieal  studies  of  travel-time  residual  strueture  as  a  funetion  of  epieentral  distanee  (e.g.  Myers  and  Sehultz,  2000). 
Our  approaeh  goes  beyond  sueh  studies  by  predieting  details  of  the  nonstationary  eorrelations,  whieh  are  diffieult  to 
resolve  with  empirieal  data  but  whieh  are  important  to  a  seismie  event  loeator.  A  erueial  test  of  our  approaeh  will  be 
loeation-validation  experiments  that  quantify  the  loeation  improvement  one  ean  aehieve  by  using  realistie 
travel-time  eovarianee  models. 
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